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We develop a general formalism to treat, in general relativity, the linear oscillations of a two- 
fluid star about static (non-rotating) configurations. Such a formalism is intended for neutron stars, 
whose matter content can be described, as a first approximation, by a two-fluid model: one fluid is 
the neutron superfluid, which is believed to exist in the core and inner crust of mature neutron stars; 
the other fluid is a conglomerate of all other constituents (crust nuclei, protons, electrons, etc.). We 
obtain a system of equations which govern the perturbations both of the metric and of the matter 
variables, whatever the equation of state for the two fluids. As a first application, we consider the 
simplified case of two non-interacting fiuids, each with a polytropic equation of state. We compute 
numerically the quasi-normal modes (i.e. oscillations with purely outgoing gravitational radiation) 
of the corresponding system. When the adiabatic indices of the two fluids are different, we observe 
a splitting for each frequency of the analogous single fluid spectrum. The analysis also substantiates 
the claim that w-modes are largely due to spacetime oscillations. 

I. INTRODUCTION 

A formalism is presented here for describing the equilibrium configurations and quasi-normal modes of general 
relativistic neutron stars taking into account superfluidity. There are both theoretical and observational reasons for 
believing that neutron stars have superfluid interiors. On the theoretical side, nuclear physics calculations |l| indicate 
that the densities in neutron stars are favorable to the formation of ^Sq neutron condensates in the neutron star 
crusts, with neutron condensates forming, as well as ^Sq proton condensates (and perhaps even kaon and pion 
condensates), in the inner core. On the observational side, there is the well-established glitch phenomenon, the best 
description of which is based on superfluidity and quantized vortices j^] . Cooling rates for neutron stars are also best 
described by superfluidity [||. 

There are several distinct ways in which a superfluid can differ from an ordinary perfect fluid. The most striking 
is that (pure) superfluids are completely free of viscous effects and are locally irrotational. The latter property is 
however compensated by the possibility for the superfluid to be threaded by quantized vortices, so that it can on 
macroscopic scales mimic the rotational behaviour of an ordinary fluid. 

In the context of neutron stars, one can distinguish three main fluid constituents: the neutrons, which can either 
belong to nuclei in the crust or be free and superfluid in the inner crust and the core; the protons, belonging only in 
nuclei in the crust and in a superconducting state in the core; and finally the electrons which behave everywhere in the 
star as a normal fiuid. As mentioned before, because of rotation, the neutron superfluid will contain vortices, whereas, 
because of the neutron star magnetic fleld, the proton superconductor is believed to contain magnetic fluxoids (at least 
if the core protons behave as a type II superconductor, which is the favored scenario |Q). Finally, the whole picture 
is complicated by the entrainment effect (misleadingly called "drag effect" in the literature) between superconducting 
protons and superfluid neutrons, by which the momentum of one constitutent carries some mass current of the other 
constituent along with it (see e.g. [|j). The theoretical description of such mixtures, including the average effect of 
vortices (and fluxoids) was given by Lindblom and MendcU j^,^ in the Newtonian context, and more recently, by 
Carter and Langlois |^,^ in the general relativistic context. 

In many astrophysical applications, the above picture can fortunately be simplified. Indeed, it has been shown that 
although the protons are superconducting, they are coupled to the "normal" electron fluid on a very short timescale 
. This implies, in practice, that the neutron star interior can be described, as a first step beyond the perfect fiuid 
approximation, as a two-fiuid model with a neutron superfiuid component and a "normal" component containing 
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everything else, i.e. the crust nuclei, the superconducting protons and the electrons. Such a two-fluid description 
was presented by Langlois et al ||lo|] in a general relativistic framework. The corresponding formalism was based on 
a formally analogous but physically different model due to Carter (see, for instance, and (ij]). The purpose of 
Carter's model was to extend to the general relativistic regime the non-relativistic two-fluid model of Landau (see 
Ip^), in which the "normal" fluid is constituted of excitations (phonons, rotons, etc.) of the condensate due to non- 
zero temperature. In the case of neutron stars, the internal temperature is believed to decrease rather quickly after 
their birth to temperatures well below the nuclear temperature scale (1 MeV), so that it will be justifled, at least for 
the global structure of the star, to ignore the temperature effects. 

The effect of superfluidity on the oscillations of neutron stars has been investigated by Lindblom and Mendell p^ ] 
in a Newtonian formalism. Although analytical solutions for simple models revealed the existence of a new class of 
modes, which could be named superfluid modes, their numerical investigation p| led only to ordinary modes almost 
indistinguishable from those obtained by a single fluid treatment, with no sign of superfluid modes. However, a 
subsequent numerical treatment by Lee p5| ], still in the Newtonian regime, was successful in extracting the superfluid 
modes. 

The main objective of the present work is, in the same spirit as |ll^ , |l5| , to investigate the influence of superfluidity on 
the oscillations of neutron stars in the general relativistic context. The description of oscillations in general relativity 
is complicated essentially by two factors: the flrst is that gravity is described by several components of the metric 
instead of the single gravitational potential of Newtonian gravity; the second is that oscillation of the star implies 
emission of gravitational radiation, and therefore the oscillating solutions arc damped (they are called quasi-normal 
modes instead of normal modes). 

Due to the difficulties inherent to general relativity, we shall restrict our analysis to idealized models. First, we 
consider only non-rotating neutron stars. Second, for the numerical application, we consider an oversimplified situation 
in which the two fluids exist throughout the star (i.e. without taking into account the crust-core boundary, neutron 
drip) and are characterized by (independent) polytropic equations of state. One of the fluids will be composed of 
neutrons, and the other will be a conglomerate made of protons and electrons, which we will call, loosely, "proton" 
fluid for simplicity in the rest of this paper. 

The changes to the dynamics of the neutron star due to the existence of two fluids will be extracted to linear order 
via an analysis of the quasi- normal modes. The linearized metric and matter variables will be decomposed into their 
respective "even" and "odd" parity components. This will generalize to the two-fluid case much of the previous work 
(e.g. [l^-|2^]) that has been done for the one-fluid case. The present work can also be viewed as a first step in carrying 
over to the general relativistic regime the pioneering calculations of Lindblom and Mendell for superfiuid 

neutron stars in the non-relativistic regime. 

In Sec . O we describe the superfluid formalism. In Sec. Ill we show how to construct equilibrium configurations. 
In Sec. |[Vf we give the linearized field equations for purely radial and non-radial oscillations inside the star and 
for the vacuum outside the star. The notation follows as closely as possible that used by |l^Jl^. In Sec. ^ the 
boundary conditions, numerical techniques, and quasi-normal mode (QNM) extraction are discussed. The techniques 
used (which are much like those of ijl^^) allow only the f-, p-, and w-modes to be reliably extracted. In Sec. ^ 
the formalism is applied to the case of two relativistic polytropes. It will be seen that the even-parity mode spectrum 
is different from that of neutrons alone, but only if the adiabatic index for the neutrons is different from that of 
the protons. In the flrst appendix we give the initial conditions that are used in the radial integration of the field 
equations. The second appendix has a brief discussion of the analytic solution used to describe the vacuum spacetime 
exterior to the star. "MTW" conventions |25| apply throughout and the units are such that G = c = 1. 



II. THE TWO-FLUID FORMALISM 



The central quantity of Carter's superfiuid formalism is the so-called "master" function. There are several choices 
for the master function. Here it is taken to be the total thermodynamic energy density —A, which is a function of 
the three scalars = —UpuP, = —pppP, and = —ppuP that are formed from , the conserved neutron number 
density current, and p^, the conserved proton number density current. The master function K{n'^ , p^ , x^) encodes all 
information about the local thermodynamic state of the fiuid and can also serve as a Lagrangian in an action principle 
for deriving the superfiuid field equations (see, for instance, [|ll],^). 

A general variation (that keeps the spacetime metric fixed) of A(n^,p^, x"^) with respect to the independent vectors 
and p^ takes the form 

5K = fipSnP + XpSpP , (1) 

where 
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= Buf^ + Apfj, , Xt^ = Cpfj, + An^ , (2) 

and 

The covectors /i^ and x^i are dynamically, and thermodynamically, conjugate to and p'^, and their magnitudes are, 
respectively, the chemical potentials of the neutrons and the (conglomerate) protons. 

The stress-energy tensor can be derived (see, for instance, Ref. [p"2|j2^ ]) from A and is found to be given by 

T,^-*<5^;+p^X. + nV. , (4) 
where ^ is the generalized pressure and is given by 

* = A-nVp-p^Xp • (5) 
The equations of motion consist of two conservation equations, 

V^n^^O , V^p^^O, (6) 
and two Euler type equations, which can be conveniently written in the compact form 

"''V[^M.]=0 , P^V[^x.]=0. (7) 

When all four are satisfied then it is automatically true that V^T^ — 0. 

Another property of this set of equations is the existence of two Kelvin theorems. Defining the (antisymmetric) 
two-forms w^i, — V[^/ii,] and Vl^i, = V[p,Xj/]j then the Lie derivative of along can be shown to be zero, and 
likewise for the Lie derivative of f2^y along . Thus, if w^^i, (H^j^) vanishes at some point on an integral curve of 
(p^) then it must do so at all other points on the curve. 

As mentioned in the introduction, this same system of equations can describe the general relativistic analog of the 
Landau model for the non- relativistic superfluid | [T^ . In this case the vector still represents a conserved particle 
number density current of the matter (which does not have to be neutrons) and /i^ is still its conjugate covector. The 
proton number density current p^, however, gets replaced with s'', which is the conserved entropy density current. 
Likewise, is replaced by 0^, the magnitude of which represents the local temperature of the fluid. 

Since we will be mainly interested in the rest of this paper by linear variations of the matter and geometrical 
variables, let us introduce now the notation we will adopt in this respect. Following |p7| , we will write the variations 
of the momentum covectors /x^ and due to a generic variation of and p^ and of the metric, in the following 
form, 

blip = A^Sp, + Bp'bn, + {bgA) pp + [SgB] Up , (8) 
5xp = Cp^Spa + A^pSn^ + (SgC) Pp + (SgA) Up , (9) 



with 



A -A .oM „M 

^dB dA dA 



Cp. = Cgp, - 2j^p^p, - 4— - jr-^npn, , (10) 



dC , dA dA 

and the terms 5gA, SgB and 5gC coming from the variation of the metric itself: 



6gA^ 



dA ^ , dA ^ , dA ^ 
on'^ dp'^ dx'^ 



(11) 



{5gB and 5gC being given by analogous formulas, where A is replaced by B and C respectively). In order to solve the 
equations for the perturbations, it will be necessary to evaluate the coefficients Apv^ Bp_i, and Cp^i, for the equilibrium 
configuration. 
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III. THE EQUILIBRIUM CONFIGURATIONS 



The background is spherically symmetric and static, so the metric can thus be written in the Schwarzschild form 

ds2 = -e'^f^'Mi^ + e-^(''Mr2 + (d^^ ^_ sin^M^^) . (12) 
The non-trivial components of the Einstein tensor are 




It is well known [ p5[ that the third of these components is not independent of the first two because of the Bianchi 
identities. 

The description of the matter is much simplified here because of the time and spherical symmetry of the underlying 
spacetime; in particular, the two conserved currents n'^ and must be parallel with the timelike Killing vector 
= (1,0,0,0), i.e. that the independent vectors are of the form 

p'" = pu^' , = nu^" , (16) 

where — t'^/|t|. Similarly, the dependent covectors are of the form 

= x{n,p)Uf, , fif, = fJ,{n,p)Uf, . (17) 

It should also be noted that x'^ = np. The stress-energy tensor also simplifies, with the non-zero components being 

r/ = A(n,p) , T::^T^ = T;^ = ^{n,p) = A + fin + xp. (18) 

The symmetry also implies some constants of the motion. Indeed, from the vanishing of the Lie derivative of /i^ one 
finds 

2i^V[^/iH + V,(i^Ai^) = . (19) 
The first term on the left-hand-side vanishes because of the equations of motion and one is left with 

- floe = t^fJ-fi ^ fJ-t ^ const . (20) 
The same procedure applies to X/i and yields 

- Xoo = f'Xf^ = Xt = const . (21) 
From Eqs. (|), (0), (|ol), and (|l|) one then finds 

fj.{n,p) ^ fiooe^"^^ = Bn + Ap , xi'^^p) = Xoce'"^'^ = An + Cp . (22) 

It is simple to verify that n^, p^, fif^ and Xti are now such that they automatically satisfy the Euler equations. The two 
conservation equations for and p^^ imply only that n and p are time- independent. Taking this with the spherical 
symmetry thus means n and p arc functions of r only. Thus, the fiuid field equations have been completely exhausted. 

On the other hand, the two independent Einstein equations — G* — 8ttT* and — 8ttT^ — still remain to be 
solved. In the case of a one-fluid neutron star, these two equations are solved simultaneously with the TOV (Tolman- 
Oppenheimer-Volkoff) equation [^ . In our case of the two-fluid neutron star, we prefer to solve simultaneously a 
system of four equations that determine n(r), p{r), A(r), and i'{r). 

The equations that govern the radial dependence of n and p are determined in the following way. Starting from 
(p^, we obtain the following equations. 
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^l' = ~^lu'/2 , ^ = -xv'li. (23) 

One can now use the generic relations (H) and (^) to express [i! and x' in terms of n' and p'. This requires us to 
evaluate the coefficients A^^y^ etc., which will be needed anyway in the following for the equations governing the linear 
pertubations. 

For the equilibrium configuration, the coefficients involving spatial indices become 

Aij — A^ij^ Bij — 13gij^ Cij — 1 (^^) 

and the mixed components, i.e. with one spatial and one time index, vanish altogether. The most complicated 
coefficients are the time-time ones. They are given explicitly by the following expressions 

.0 , ^dB ^dA n ^dA . OA 
op"^ on^ op'^ ox'^ 

„ ^dB ^ dA dA 2 
an^ ow^ ox'^ 

^o-C + 2—p^ + A—np+—n^. (25) 

In the three coefficients above, after the partial derivatives are taken, then one sets x'^ = np. 

The expressions ( p3[ ) can now be rewritten as differential equations that determine the radial profiles of n{r) and 
p(r): 

Ay + By + ^{Bii + Apy ^ , Cy + Ay + ^{An + Cpy ^0 . (26) 

The functions vir), A(r), n(r), and p{r) can now be determined by solving these two equations in conjunction with 
the two independent Einstein equations, which can be written in the form 

1 — f — e'^ 

A' = 8nre^A{n,p) , i^' = h 87rre^^'(n,p) . (27) 

r r 

This section is ended by recalling that a smooth joining of the interior spacetime to a Schwarzschild vacuum exterior 
at the surface of the star (the radial value of which we define to be R) implies two things: (i) the total mass M of 
the system is given by 

M = -An / dr A(r) (28) 



and (ii) the total radial stress must vanish at the surface of the star, which implies ^{R) = 0. In this work we will 
only look for background solutions that satisfy n{R) = and p{R) — 0. For the master function used later, these two 
conditions guarantee that ^'(i?) = and also lead to A(i?) = 0. In view of (|27|), requiring a non-singular behaviour 
at the center of the star will impose that A(0) = 0, and consequently that A'(0) and i''(0) must also vanish. This in 
turn implies, in view of (26), that p'{0) and n'(0) have to vanish as well. 



IV. THE LINEARIZED FIELD EQUATIONS 

The metric and fluid perturbations will be decomposed on the basis of spherical harmonics Yi"^{9,(p) (where I is 
the angular momentum number and m is its projection onto the z-axis) [E8[. Because the background is spherically 
symmetric, there is no real loss of generality by restricting the study to the modes m = 0, i.e. by considering 
perturbations that do not depend on the (p coordinate. This means then that the basis on which the perturbations 
are actually decomposed is simply that of the Legendre polynomials Pi{d) = Yj^^^{0, (p). 

It will be seen to be convenient to distinguish between the so-called "even" (i.e. (— I)'-parity) components and the 
"odd" (i.e. (— l)'"'"^-parity) components. For the even-parity modes a time dependence of the form 

0(r, i) = C'(r)e*'^* (29) 
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is assumed. In general, oj is complex. The same time-dependence for the odd-parity modes will not be imposed, 
because it will be shown later that there are no odd-parity pulsational modes for the fluid, but just differentially 
rotating ones (see also for a discussion of the same result for one-fluid stars). 

Purely radial oscillations correspond to Z = 0. However, the system of equations written with general I in mind 
must necessarily exclude the case of ? = 0. This makes it necessary, at the appropriate time, to consider the radial 
oscillations separately. 



A. The Linearized Metric and Matter Variables 



The form of the linearized Einstein Equations used is 



where 



(30) 



(31) 



and T = g^'^T^i, is the trace of the stress-energy-momentum tensor. The fluid equations to be linearized are the two 
conservation and two Euler equations, of which the latter two take the simple looking forms (for both parities) 



(32) 



In the decomposition of the perturbations into even and odd-parity components, it can be seen that the ten linear 
perturbations of the metric will be divided into respectively seven and three components. Using then the four gauge 
transformations of general relativity, i.e. the coordinate transformations, three acting in the "even" subspace and 
one in the "odd" subspace, we can reduce the description of the metric perturbations to four "even" components 
and two "odd" components. A convenient choice is the so-called Regge- Wheeler gauge ||2^, in which the even-parity 
components of the metric perturbations are 






























(33) 



whereas the odd-parity components are 
^3m<^ = 



hQ{r,t)sm9-^ hi{r,t)sm9-, 



ho{r,t)sme-§g 

hi{r,t)sm9-§g 





(34) 



To facilitate further manipulations of the linearized Euler equations, conservation equations, and STp^, and 
are rewritten as a product of a magnitude with a unit timelike vector so that 



Pp = PVp 



(35) 



{u^Up = — 1 and v^Vp = — 1), where contrarily to the background case Up and Vp are not aligned in general. For our 
background, the spatial components of the perturbations of the momenta (recall Eqs. (H) and ^) are given by (true 
for both parities) 



S^i = BnSui + ApSvi , Sxi — CpSvi + AnSui 



(36) 



The time components are more difficult to obtain, but a straightforward calculation, using Eqs. (|), (|) and (|ll|) and 
inserting (p3|) and (53), leads to 



Al^v + ^l^n + - (Sn + Av) -ffoPje*"* 



(37) 



and 
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Sxo 



C^Sp + A°oSn +-{An + Cp) H^Pie''^' 



(38) 



for the even-parity components and S/io = ^xo = for the odd-parity components. 

The hnearization and solving of the conservation equations, V^n^ — and V^p^ — 0, are made easier (for even- 
parity modes) by introducing S^l^ and S^p, which are the Lagrangian displacements, respectively, for the neutron and 
proton number density currents. Now, 



c\ f) d (9 



dt 



dt 



(39) 



where r„ (rp) is the proper time as measured by an observer whose four- velocity is {'^^)- A decomposition of the 
Lagrangian displacements into spherical harmonics yields 



5C = -r'-2K(r)^P;e-* 



and 



r _ ^-A/2J-1 



r'-^iyp(r)P,e^- , ^ -r'-'Vp{r)-Pie 



d 



de 



(40) 



(41) 



Insertion of the Lagrangian displacements into the conservation equations castes these equations into a form that 
can easily be integrated, the results being 



An(r) 
n{r) 



/ + 1 1 

-^Wn + -Wn' 



1(1 + 1) 1 



for the Lagrangian variation in the neutron number density, and 



Ap(r) 
pir) 



(42) 



(43) 



for the Lagrangian variation in the proton number density. These are related to the Eulerian variations 5n[r) and 
5p(r) via (see [psl, pg. 691) 



Ap = 5p + p'e^^^^r^-^Wp . 

For the odd-parity modes, Sn{r) = and 5p{r) — 0, and the only non-zero components for (5u^ and Jw^ are 

8Pi BPi 
6u^ = e-^/^Unir, t)sme—^ , Sv^ = e-^/^Upir, i)sin(?— f . 

Uu Ou 

It can be shown that the conservation equations are automatically satisfied in this case. 



(44) 



(45) 



B. Equations for Radial Oscillations 



When / = more gauge freedom exists to eliminate some of the metric perturbations. This freedom will be used to 
set Hi = and K = Q. Since the Legendre polynomial Pq is just a constant, then all of the odd-parity perturbations 
and the Lagrangian displacements (5^^ and (5^p vanish. The only remaining perturbations are TJq, H2^ Wn, and Wp. 
The remaining field equations consist of four linearized Einstein equations and two linearized Euler equations for the 
fluid. 

An independent set of equations contain two for the metric, which are 

H'a = Airre^ (p^C^ + 2npAl + n^B^ - 2* - H2- 
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S7rp>/2 

H2 ^ Sire^/^XpWp + fxnWn) , (46) 
and then two more for the fluid perturbations, 

, ,2 / „(i'-A)/2 

r 



^j2 / (iy-A)/2 

U'^/'[C',p + A'on]H2) +^XHI,. (47) 



The two equations for the metric perturbations can be used in the two for the fluid perturbations to reduce the 
total system to two coupled second-order differential equations for Wp and Wn- This generalizes to our case the 
one-fluid results of Chandrasekhar who was able to go further and demonstrate that the final equation for the 
perturbations is of the Sturm-Liouville form. Whether or not our system forms such a self-adjoint problem remains 
to be established. 



C. Equations for Non-radial Oscillations 



To obtain the equations for the even-parity oscillations, we have followed the formulation of Detweiler and Lindblom 
[p^ , the main difference in our case arising from the existence of two fluids instead of one. The procedure to obtain 
the equations is the following: the first step is to write explicitly the linearized Ricci tensor on the left-hand-side of 
( ^ in terms of the quantities introduced in (js^); the second step is to write the right-hand-side of (|3^) in terms 
of the matter variables Wp, Wn, Ki and Vp, as well as the metric variables, using (p6|)-(^^. One then obtains six 
relations. In the case of a single fluid, there are enough equations. However, in the case of two fluids, additional 
equations are necessary. They will be given by the equations of motion for each fluid. The conservation equations are 
automatically satisfied as a result of our choice of Lagrangian variables. What remain are the Euler equations. This 
will provide us with four equations, corresponding to a radial component and an angular component for each fluid. 
The whole system of equations is now redundant because of the Bianchi identities. 

After tedious calculations (performed by hand, and independently with the aid of Mathematica), we found that the 
Einstein and fluid equations yield a constraint equation 



2-1 -f 
-2e 



^(1 



„2A 



2c^ 



-A 



1 - - (1 - e"-^) - 47rr2* 



Hi 
K 



(48) 



having no radial derivatives, and the rest, which are first-order in the radial derivative, given by (using the notation 

T)0 _ no no _ /l02\ 



Hi' = -Ho 
r 



A' 



l + l 



Hi + —K - 167r— {pnVn + XPVp) 



(49) 



K'^^+^lplHi 
r 2r 



l + l 



gA/2 

K-8n [finWn + XpWp 



(50) 
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' + 1 "'it.. 
+ — W^n- 

r n 



(51) 



(52) 



with 



— 



lx„ + — 

r 2 



6^/2 l{l + 1) u;^ 



K 



_e(^--)/^—n {BnWn + ApWp) - A^e^^+''^'^^ {unW,, + XpWp) 



= -(A-'>)/2 



^!L(BO'nW,,+A','pWp) + {^ 



n 

2r r 



{B^nWn + A^opWp)] 



-X 



p 



+XP 



2 



PX (i-^^')-p' (CSp + ^SJn) 



Xp[j-Y^]-{C',P + Xn)p' 



K 



e'^'^p' {C^^pVp + AlnVn) 



„e(^-'^)/2^^p {iZpWp + y4nW„) - 47re(^+'')/2 ^^^^^ ^^^^^ 



-^-{c^,'pWp + AUWn) + r-^ 



2p' A' - 1/' , p" 



2r r 



Xn=n 



-Mi?o + e'^/^u^ {BnVn + ApVp) 



e^'^-^)l^-{BlnW^+AlpWp) , 



(53) 



(54) 



(55) 



Xri 



-xHo + e-^/W {CpVp + AnV„) 



{C°pWp + A°onWn) 



(56) 



Note that H2 has disappeared. This is because the Einstein equations (for the present case of ^ > 0) impose H2 = Hq, 
in contradistinction to the purely radial case. 

Let us emphasize that the form of our system of equations is quite analogous to that given in for the one-fluid 
case. The main difference in our system is the doubling of the number of equations corresponding to the fluid degrees 
of freedom (i.e. the W's and V's). Note also that an explicit coupling between the two fluids is manifest in the latter 
equations when the coefficients Aq or A are non-zero. If they are zero, this means the two fluids are independent but 
they remain coupled indirectly through the gravitational field, i.e. the equations governing the degrees of freedom of 
the fluids include a dependence on the metric perturbations. 
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For odd-parity perturbations the Euler equations yield 



ApUp + BnUn = 
whereas the non-trivial linearized Einstein equations are 



CpUp + AnUn = 



A' 



2r 



2r2 



hi 



hi 



1 
2e^ 



h' 



-ho 



47r(^' + A)ft.i 
j^' -X' 



2v' _ 



-ho - 



h'l 



IGtt 



v' + \' 



hi 



2 

/(/ + 1) 



hi] =Q 



ho = 



+ finUn ] +8t:{^ + A) ho 



(57) 



(58) 



Recall that odd-parity modes automatically satisfy the conservation equations. The Euler equations thus imply that 
Un and Up both vanish, which means at this order the fluid only differentially rotates, and does not pulsate. This 
is similar to what happens in the ordinary one-fluid case JTof . If pulsations are considered, then the linearized Euler 
equations imply that the fluid perturbations vanish, and then the gravitational degrees of freedom are decoupled from 
those of the fluid. 



D. The Linearized Equations Outside the Star 



Outside the star, the fluid variables are all zero. The problem is then reduced to the metric perturbations alone, 
which can be conveniently obtained from the so-called Zcrilli function Z |^^, which satisfies 

^ = {VAr)-u;^)Z, (59) 

where r* = r + 2Afln(r/2Af — 1), M is the stellar mass, n — {I — + 2)/2, and the effective potential is 

1 — 2M/r 

yz{r) = TTTT^ i'^ri'^in + ly + Gn'^Mr^ + ISnM^r + 18M^) . (60) 

As r ^ r* — > oo, then the asymptotic form of Z is 

Z ^ A^^e'^-' + Aoute-'''''' . (61) 
The relations [^,0 that determine Hq and K from Z are 

,^ _ a{r)Ho{r) + {h{r)-k{r))K{r) 



and 

dZ (r*) 



dr* 
where 

, , nr + iM 

a(r) 



^ K{r) ~ g{r)Z {t*) , (63) 



h{r) = 



^2^2 _ (j^j^ l)M/r ' 

nr{r - 2M) - u'^r'^ + M{r - 3M) 
(r - 2M) {uj^r^ - {n + l)M/r) 
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n{n + l)r^ + 3nMr + 6M^ 
r2(nr + 3M) 

-nr^ + 3nMr + 3A'P 
(r-2M)(nr + 3M) ' 

^2 



fc(r) = - 



r-2M 



(64) 



(It should be mentioned here that b{r) suffers from a typographical error in p?!, but was subsequently corrected in 
[^6|.) The remaining metric perturbation Hi can be obtained from the vacuum form of the constraint equation (|4^), 
i.e. 



3M + ^{l + 2){l - l)r 



Ho = - 



2 



-(? + 2)(l - l)r - wV^e"'' - -e^M{3M - r) 



(65) 



E. The role of vortices 



A superfluid differs from an ordinary fluid in that it is locally irrotational. It can nevertheless closely mimic rotation 
via the formation of an array of quantized vortices that threads the fluid. Such an array has not been explicity included 
here. Since our initial configuration is not rotating, it contains no vortices. When there are no vortices, the superfluid 
is said to be in a Landau state ||3l[| . A further restriction must be imposed to obtain a strict Landau state. This 
restriction is that the vorticity — V[^/i^] vanishes. This is already automatically satisfied by the background 
configuration. 

If one wished to impose the irrotationality condition on the perturbations as well, then it can be shown that this 
would restrict the set of allowed perturbations by imposing the constraint Hi (r) = in the case of even-parity and 
BnUn + ApUp in the case of odd-parity perturbations. However, a small perturbation at the scale of the whole 
star will be felt at a mesoscopic level as a huge perturbation in the sense that many vortices will form spontaneously. 
This can be illustrated for example by a rough estimate of the critical angular velocity for the formation of one vortex. 
For a Newtonian laboratory superfiuid, this value is given by (see e.g. ||3l[] ) 

where a is the size of the vortex core. For a neutron star, a ^ 10^^^ cm and R ^ 10 km, so that fici ^ 10"^'* rad/s, 
which is extremely small. This means that, even for small perturbations of the star, one expects a huge number of 
vortices so that the superfiuid will mimic an ordinary perfect fiuid on macroscopic scales. It thus makes perfect sense 
to treat the superfiuid as we have done here. 



V. BOUNDARY CONDITIONS, NUMERICAL TECHNIQUES, AND QNM EXTRACTION 



A. Boundary Conditions at r = and r = R 



In Sec. VI numerical solutions to the even-parity system of equations ([48[)-(p4D are constructed. An essential 
element of this construction is to know the behaviour of the fields near r = 0, because of a general problem of 
spherical coordinates that one cannot start a numerical integration precisely at r = 0, but only from a nearby point. 
As outlined by the values of the fields near r = can be obtained by expanding each field in a Taylor series, and 
then using the field equations to determine the coefficients of these expansions. The results, which are presented in 
the first appendix, show that one need only specify the set of values {if (0), W„(0), Wp(0)}, for instance, and then the 
remainder, i?i(0), Xji(O), and A^p(O), are determined. All of the second derivatives, iJ"(0), K"{Q), etc, are likewise 
determined by {X(0), W„(0), Wp(0)}. 
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At the surface, the story is somewhat different. In the case of a single fluid star, the boundary condition is that the 
pressure should vanish at the physical surface of the star. Since the fluid has moved with respect to its equilibrium 
value, the physical surface is displaced with respect to the background star surface and the appropriate condition is 
thus that the Lagrangian variation of the pressure vanishes at the surface of the star. The same reasoning cannot 
be applied here because we have two Lagrangian variations to consider, one along the neutrons and another along 
the protons. However, in our particular example where the two fluids are independent, one can generalize easily the 
single fluid result because the pressure is separable, in the sense that it can be decomposed unambiguously into two 
pressures corresponding respectively to each fluid, 4" = 5'„ + 5'p (with = A„ + /in, etc.). The physically meaningful 
boundary condition is then to impose that the Lagrangian variation of each pressure should vanish at the star surface. 
Let us express this condition in terms of our variables. The Lagrangian variation of for instance will be given by 

A*„ = nAi.1 = nB^An . (67) 

Substituting the expression (|2|) and using H2 = Hq as well as (^ij) , one gets the simple relation (like in the one- fluid 
case of course) 

A*„ = -e-^/V'X„ , (68) 
and the corresponding result for A'^p. Our boundary conditions will thus be X„(i?) = and Xp{R) = 0. 



B. Numerical integration for quasi-normal modes 

One finds in the numerical integration the closestparallels between the formalism presented here and what has been 
done for the perfect fluid with one constituent ||l6|,|l3 . Given a background configuration, the three main components 
of the numerical procedure are (i) integrate inside the star, (ii) integrate outside the star, and (iii) match the interior 
and exterior solutions at the surface of the star, and find those lu (corresponding to quasi-normal modes) for which 
there are only outgoing gravitational waves at spatial infinity. We use the dgear subroutine available from the IMSL 
library to numerically integrate the background equations, and the linearized equations inside and outside of the star. 

The integration inside the star has itself two separate parts. The perturbation equations (|4^)-(|5^) in the interior 
can be written as the matrix equation 

dY 

— = Q.Y, (69) 

where 

Y = {Hi,K,Wn,Wp,Xn,Xp} (70) 

is an abstract six-dimensional vector field. The 6x6 matrix Q depends on Z, lu, the background fields, and the 
various coefficients A, Aq, etc. As mentioned earlier, the results presented in the first appendix indicate that only 
the set {K{0),Wn{0),Wp{0)} needs to be specified and then the remaining quantities i?i(0), A'„(0), and Xp{0), and 
all of the second derivatives at r = 0, are determined. At the surface of the star, there are only the conditions that 
X„(i?) = Xp{R) = 0, so that one must specify the set of values {Hi{R),K{R), W„(i?), Wp{R)}. 

For a given background, the first part to the interior integration consists of choosing three arbitrary values of 
{K(0), Wn{0), Wp{0)}, and then starting the integration of Eq. ( |69| ) at small rg, where 

Y(ro)=Y(0) + V(0)r2 . (71) 

(Both vectors Y(0) and Y"(0) can be constructed from the results presented in the first appendix.) The integration is 
terminated at r = R/2. This process is repeated two more times, each with a different choice for {^(0), Wn{0), Wp{0)}, 
to get three linearly independent solutions Yi, Y2, and Y3, and hence a general solution in the domain Q < r < R/2 
of the form 

3 

Y(r)=^c.Y,(r), (72) 
1=1 

where the Ci {i — 1, 2, 3) are constants to be determined. 
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The second part of the interior integration starts from the surface of the star. There are four Hnearly in- 
dependent solutions Y(r) that can be constructed. These are built by choosing four arbitrary sets of values 
{Hi (R), K{R), Wn{R), Wp{R)} and then integrating backward from Rto R/2 four times to generate a general solution 
in the domain R/2 < r < R oi the form 

7 

Y(r)=^c.Y,(r) , (73) 

where the c; (i — 4, 5, 6, 7) are constants to be determined. 
At r = R/2, the two general solutions must be equal: 

3 7 

^c,Y,(i?/2) = ^c.Y,(i?/2) . (74) 

i=l i=4 

Each value of Yi(R/2) (i = 1, 7) is of course known. Hence, picking arbitrarily one of the q and giving it some 
value, then Eq. ([74|) represents six equations that can be used to determine the remaining six c;. Having obtained 
the Ci, Y(r) is thus known throughout the star. 

Next, the metric perturbations outside the star must be determined. A priori, since all perturbations are known 
throughout the star, the boundary values Hq(R) and K{R) obtained from the interior solutions are sufficient as initial 
data to determine all the metric perturbations outside the star, embodied in the particular Zcrilli function Z^eg which 
satisfies, according to (|63) and (|63|), 



dr* 



^reg 



_ aiR)Ho{R) + (6(i?) - fc(i?))X(i?) 
~ h{R) - k{R)g{R) ' 

= K{R) - giR)Zreg . (75) 



However such a procedure will give in general a mixture of ingoing and outgoing gravitational radiation, whereas we 
wish to determine the specific values of u for which there is only outgoing gravitational radiation. For the conventions 
used here, this implies that = in the asymptotic form of the Zerilli function (cf. Eq. (|6^)). It has long been 
known, however, that this is problematic numerically, because stable stars undergoing quasi-normal mode oscillations 
must have (for the conventions used here) Im(aj) > and this implies that the exponential multiplying Aout in Eq. 
( |6l| ) diverges as r* — > oo. We use two different techniques to circumvent this difficulty. One will enable an accurate 
determination of the so-called w-modes, and the other an approximation of the f- a nd p -modes. 

To determine the w-modes we use the Leaver series [^2[^ in a manner like that of |34-36|. It is, albeit indirectly, an 



analytic representation for Zout, the Zerilli function that describes an outgoing wave at spatial infinity. (The details 
of this series are given in the second appendix.) Like Zreg, Zout depends only on I, lo and M , so that for given I and 
M the only free parameter is w. The special values of w that correspond to quasi-normal modes are extracted by 
solving for the (in general, complex) roots of /(w) = 0, where 



\ I ^Z^^q 

f[UJ) = -— ^ 

Zreq dr* 



1 dZ„ 



dr* 



(76) 



We use MuUer's method |Q to determine the roots. 

In principal, the procedure just outlined should succeed in determining all the quasi-normal modes. In practice, 
it only works well for the w-modes. The difficulty lies in the fact that the ratio Im(a;)/Re(a;) of the fluid modes 
are typically so small (< 10~^; see also [^) that the iteration of MuUer's method is not possible. (In practice, 
one can only obtain the f-modes by using very accurate initial guesses.) To determine the f- and p-modes a direct 
integration is done, with Zreg and dZreg/dr* (cf. Eq. ([75[)) as initial conditions, to obtain the Zerilli function. 
The problem of the divergence of Z is overcome by looking only for resonant solutions that have real uj. The true 
quasi-normal modes are not produced, but accurate numbers for the real part of the frequencies should be obtained 
since Im(w)/Re(tj) < 10"''. The radial profiles of all the fields should thus also be representative of the true quasi- 
normal modes. In general, direct integration will yield a Z that corresponds to both incoming and outgoing waves at 
spatial infinity. The resonant frequencies that correspond to an outgoing wave alone are extracted via a "graphical" 
technique that maps out log \Ain\ vs iv. Those values of co that cause deep minima to occur in log \Ain\ correspond to 
the resonant frequencies. 
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C. Convergence Tests and Accuracy 



Convergence tests have been performed to determine the rehabiUty and accuracy of the numerical results. The 
tests discussed here are for the model three presented in Table I, which is for a star with neutrons and (conglomerate) 
protons, each fluid behaving as a relativistic polytrope (cf. Eq. ([79|)). Also, we have considered only I = 2 modes. 

The accuracy of the results for the w-modes is determined by calculating the changes in the w-mode frequencies for 
different step sizes. The integrations have been performed using (scaled) step sizes of Ar = 10~^ (which corresponds 
to Ar/R ~ 10""*) and Ar = 5 x 10~^. The fractional change in the real part of the frequency is defined by 

Re(Ac^) ^ Rc(cj2) -Re(^i) 

Re(c^) ~ Re{uji) ' ^ ' 

where cui is the frequency calculated using the first step size and 0J2 is for the second. A similar definition is employed 
for the imaginary part of u. The results for the first six w-modes are summarized in Table II. 

Another test has been done by changing the matching point inside the star (cf. the numerical integration discussion 
of the previous subsection). Two different matching points, R/2 and i?/4, have been used. It is found that the 
fractional change in the frequencies is Aw/uj ^ 10~* — 10~^, which is comparable to the changes of Table II. One 
can be confident, then, that the w-mode frequencies produced by our numerical scheme have at least three significant 
figures accuracy. 

The accuracy of the f- and p-modes is determined by varying the truncation point rtmnc which is the final grid 
point of the direct integration that determines the Zerilli function. As discussed earlier, the resonant frequencies are 
located by finding the deep minima of the incoming wave amplitude Ain, which is obtained from 



luj dr* 



(78) 



The tests reveal that good convergence is achieved when rtmnc > 50i?. It was also found that the resonant frequency 
corresponding to the first deep minimum can be used as a good initial guess in the Leaver's series method, and agrees 
very well with the real part of the f-mode obtained from the Leaver's series {ujM — 0.130 -I- 6.52 x 10~^i). A value of 
Ttrunc — 300R wiU bc used in the next section. 



VI. APPLICATION: TWO RELATIVISTIC POLYTROPES 



We will now apply our formalism to the concrete example where each fluid behaves as a relativistic polytrope. Let 
us take a master function of the form 

A(ri^,p^) = -m„n - (T„n'^" - m„p - app'^" , (79) 

where we have assumed for simplicity the same mass m„ (i.e. the neutron mass) for both fluid particles. This master 
function is not only independent of x'^ it is also separable. An immediate consequence is that A — A!q = 0. The other 
relevant coefficients are given by 

SS = cT„/3„(/3„-l)n^"-2^ (80) 

Co"-fTp/3p(/3p-l)/''-^ (81) 

The two fluids are also assumed to be in chemical equilibrium for the background, which means that the condition 
/i = X is imposed, and therefore n and p are linked by an algebraic relation valid throughout the star. Thus, only one 
central particle density needs to be specified in order to determine the star model. It should be noted that chemical 
equilibrium is entirely consistent with Eqs. and (|2^) (or Eq. (^)), since these equations themselves imply that 
if /i = X one point in the star then /i = x at all points. 

The discussion will be limited to two particular choices for the free parameters (t„, CTp, /3„ and which are the 
two different sets listed for models one through four in Table I. In Fig. 1 we show the dependence of the mass M on 
the total central density, i.e. Uc = uq +po, using the "cr" and values of models one through four in Table I. Fig. 
2 contains the radial profiles of n and p for models one (left graph) and two (right graph). One can see in Fig. 1 the 
existence of maxima at around ric = 2.5 fm~^ for both sets of "ct" and values. The branches of these curves to 
the left of the maxima represent those configurations that are stable to radial perturbations. Thus, the two sets of 
values for no and po chosen for models one and two (to be used in the quasi-normal mode analysis below) correspond 
to stable configurations; models three and four are unstable. 
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A. The w-modes 



The w-modes were first discovered by Kokkotas and Schutz p3q | for one-fluid stars. They are due primarily to the 
oscihations of spacetime itself, coupling only weakly to the fluid of the star. (See, for instance, Andersson et al 
who extracted w-modes using an Inverse Cowling Approximation, wherein all the fluid motion is frozen out.) Another 
characteristic property of w-modes is that they are strongly damped (Im(a;) ^ Re(a;)). Lein et al have found 
a branch of strongly damped w-modes, which are denoted wn, that are similar to the quasi- normal modes of black 
holes. 

In Fig. 3 we plot, for I = 2, several w-mode frequencies for models one and two, and in Table III we list the explicit 
values for the real and imaginary parts for the first six modes of Fig. 3 for both models. The wn-modes are the two 
whose frequencies satisfy Im(a;) > Re(u;). There are no qualitative differences for the w-modes between our models 
one and two, nor are there such differences between both models and what one finds for a one-fluid system. Figs. 4 
and 5 are radial profiles (for the wi- through W4-modes listed in Table III) of Wn (which also describe Wp since it 
is indistinguishable from Wn for these cases) for models one and two, respectively. Again, there are no qualitative 
differences between the two models for W„ (or Wp). 

It is also important to point out that Figs. 4 and 5 demonstrate that the neutrons and the protons move in "lock- 
step" with each other for each w-mode frequency and for both models (the density variations show the same). This 
can be considered to be further evidence for the idea that w-modes are largely spacetime oscillations, because below it 
will be seen that the obviously fluid oscillation f- and p-modes allow for counter-moving motion between the neutrons 
and the protons. 

B. The f- and p-modes 

The f-mode, or fundamental mode, may be regarded as the "simplest" non-radial mode. It has no radial nodes and 
it usually reaches a maximum value near the surface of the star |4l| ]. Physically one may understand this mode as a 
surface wave that maintains the total volume of the system. The p-modes, or pressure-modes, have as their restoring 
force pressure and they behave like sound waves. Typically there is an infinite number of p-mode frequencies. The 
radial profile of the first p-mode, which will be denoted pi, will contain one radial node; the second mode, denoted 
by p2, will have two; and so on. Unlike the w-modes, there are qualitative differences between models one and two 
for the f- and p-modes. 

Fig. 6 is a graph of log|v4i„| vs Re(a;M), with I = 2, for model one and also for a star composed solely of neutrons 
behaving as a relativistic polytrope (with the same neutron parameter values as model one). Reading it from left to 
right, the first deep minimum corresponds to the f-mode, with the next corresponding to the first p-mode, and so on. 
One should notice that the deep minima for the one-fluid star are virtually indistinguishable from the star with two 
fluids. 

Fig. 7 is a similar graph, but for model two of Table I. The obvious difference with Fig. 6 is the doubling of the 
f-mode, and (near) doubling of the p-modes. For a more precise numerical comparison. Table IV lists the real parts of 
the f- and pi-mode frequencies for model one, and their "doubled" counterparts from model two, which are denoted 
fo ("o" for ordinary), fg ("s" for superfluid), pio, and pig. 

The new, i.e. superfluid, modes represent the case where the neutrons and protons no longer move in "lock-step" 
with each other. This is illustrated most dramatically, perhaps, by Figs. 8, 9, and 10. Fig. 8 contains the radial 
proflles of An/n (or Ap/p, since they are identical in this case) for the {I = 2) f- and pi-modes of model one; Fig. 9 
gives both An/n and Ap/p for the {I = 2) fo- and fg-modes of model two; and. Fig. 10 gives An/n and Ap/p for the 
{I = 2) pio- and pig-modes of model two. For the ordinary modes the two density variations increase and decrease on 
the same intervals, but for each of the superfluid modes the two density variations are clearly out of phase, i.e. the 
neutron number density increases when the protron density decreases, and vice versa. 

Figs. 11, 12, and 13 contain (for I = 2) radial proflles of Wn and Wp for models one and two. Fig. 11 is for model 
one, and not surprising there is no difference between the radial flow of the neutrons and the protons. However, Figs. 
12 and 13 are for model two, and clearly evident is counter-motion between the neutrons and the protons for the 
superfluid modes, i.e. the neutrons are flowing in when the protons are flowing out, and vice versa. 

Finally, with Figs. 14-17, one can contrast the effects of the superffuid with the ordinary fluid modes on the 
spacetime perturbations Hi and K. We have also considered higher I values, and have found superfluid modes. Table 
IV lists the higher I frequencies for the f- and pi-modes for both models one and two. 
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C. W-, f-, and p-modes for Radially Unstable Configurations 

The W-, f-, and p-modes have been investigated for configurations that are unstable to purely radial oscillations, 
that is, they are located to the right of the maxima of Fig. 1. These are models three and four of Table I. The general 
features of the quasi-normal modes discussed above for models one and two carry over to models three and four. Thus 
even stars on the unstable branch contain the counter-moving modes. 

VII. CONCLUDING REMARKS 

The main purpose of this work was to build a formalism appropriate to study the quasi-normal modes of general 
relativistic neutron stars with superfluidity. We have described the matter content of a neutron star in terms of 
simply two components, one corresponding to a neutron superfluid and another corresponding to a proton-electron 
conglomerate fluid. This approximation makes sense since protons and electrons can be seen as locked together, from 
their electromagnetic interaction, on bulk motion scales so that they move in "lock-step" with each other as the star 
oscillates. We have ignored effects due to the crust, as well as electromagnetic effects that would lead to, for instance, 
the magnetic fields that neutron stars are known to possess. 

We have performed a numerical investigation based on the system of equations that we obtained for linear pertur- 
bations both of the fluids and of the metric. We have found numerically general relativistic superfluid modes, i.e. a 
new set of modes characterizing the dual nature of the star matter content which doubles the number of degrees of 
freedom in the matter. In the Newtonian regime, these modes had been revealed analytically for a simplified model 
by Lindblom and Mendell [p^ , p4| j5||^ and found numerically by Lee p^ . 

A distinguishing characteristic of these new modes is that they are counter-moving. Also, when the adiabatic indices 
of the two constituents are equal then there is a degeneracy between these new modes and their analogs in one-fluid 
stars. Lindblom and Mendell argued that the superfluid mode frequencies should be larger than the ordinary mode 
frequencies. Our work verifies this conclusion for the general relativistic regime: Table IV clearly shows that fo < fg 
and pio < Pis- We have also verified the inequalities for the higher node p-modes. 

Finally, our results also tend to support the claim that w-modes are due mostly to spacetime oscillations, since 
there is very little qualitative difference between our models one and two. The obviously fluid f- and p-modes, on the 
other hand, were clearly distinguished between models one and two, exhibiting the counter-moving modes. 
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VIII. APPENDIX I: RADIAL INTEGRATION INITIAL CONDITIONS 

The background quantities will be expanded near r = in the form 

q{r)=q^ + \q2r^ + \qir^^O{r^). (82) 

The field equations imply that the first and third derivatives of each field must vanish at r = 0. We shall use an 
overbar and an overhat to designate the zeroth order and the second order respectively, so that, for instance, 

n = no , n= , W = fiQUo , \ (^^""2 + ^^2no) ■ (83) 

By comparison with a Taylor series expansion, one recognizes that 
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qo = q{0), q2 = q"{0), q^ = lq('-\o). 

D 



(84) 



Now, the lowest order of the perturbation equations yield the following three constraints (where it is useful to know 
that Hq{0) = K{0), Vn{0) = -Wn{0)/1, and Vp{0) = -Wp{0)/iy. 



^n(O) = -^WK{Q) - W'^'^n^Bln + —e-'^I^BnA W„(0) 



(85) 



^p(O) = -^XPm) - [e'">/^P2C°oP+—e-''°^^Cp^j %(0) 
- (|e'^°/2p2^+ ye-''°/234^^ W„(0) , 

For the next order, it is useful to recognize that the background Einstein equations imply 

A"(0) = -^Ao, ^"(0) = y (34-0 - Ao) . 



It is also useful to know that 



and 



where 



h;;{o) = k"{o) + Qo 



1 + 3 



1 + 3 



K"(0) = Q„ - ^^^VF„"(0) , V/{0) = Qp-j^^^W/{0), 



l{l + l) 



{l + 2){l-l) 
_ /27rZ(/ + l) 

"V 3 



87re-^°/2 (X„(0) +Xp(0)) - (^-^A^ + u;^e-'°^ K{0) 
(3*0 - Ao) - w^e-^") fl'i(0)l , 



(86) 
(87) 

(88) 
(89) 
(90) 



(91) 



'-0 



[e-""'/%(0) +P2 (^Wn{0) +C^WpiO) 



(92) 



Qp — 



_Bl_ 



( -^8 



e-'^''/2Xp(0) +P2 (CoVWp(O) + ^gnW„(0)) 
Xn(0) + n2 (^W^p(O) + SO^W„(0)) 



-yo/2 ■ 



The next order of the perturbation equations then yield the following: 



1 



2 -xJ'{o)---fmK"io) + - 



Wn"iO) 



(93) 
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i{i + iy 



2 
2 



Lu I Bn^Qn + ApnQp 



714 



47rAo 



g-fo/2 I I 

i/2X„(0) + -lmK{0) + -JmQo 



Apn 



Wp{0) , 



ni + 



47rAo 



3 

"2 ) Alp + H'iAlp + 



Ul 



{^Anp— 



p-vojl 1 1 

^X/'(0) - ixP^"(0) + \ 



+ 



2 

z^"(0) 
2 

z."(0) 



;(z + i) 



Anp 



47rAo \ -Q- 



^ i^2^p(0) + ixP^(0) + ixpQo 



P4 



Wp(0) - 



Pi + 



47rA( 



3 

P2 ] A^n + p2A^n + ^ (Anp 



1677 „ . 1 



[/mW„(0) + xpWp{0)] - Stt (7InQ„ + xpQp) + -Qo , 



1 + 2 



8 

r / _|_ 2 1 

W„"(0) + e""/^ -^n2Aln + 2TTjmxp + -e-'-'u^Apn 



■7InJfi"(0)+e"»/2 



/ _(_ 2 1 

—^n2B^n + 2Tr{jmf + -e-'^'co^Bn^ 



\IJ,n 4 / 



1 3 



W^,"(0) 



+ |e''o/2 /^^/„2 _ (/ + 2)n4 - in2Z^2 - ^ Aon2 ) Bgn - Zn2Sgn - 'iTr]mim+ 
^ (47r7Zn)' Ao - non2 (SS)^] + a;2e-'^°/2 + 1^.2 + yAo) - B^?] 1 1 



,'^o/2 



(jjjfl \ 4:77 \ "' 

—^712 - + 2)714 - 2"'2'^2 - -^Aon2 j A^p - ln2Alp - 4TTjmXP+ 



^ + ^1^2 + yAo^ Am -Apn^ Wp{0) , 



Xp"{Q) 



1 + 2 



;a+l)„.„/2_ 



'^°/2^iJi"(0) + e''''/" 



7 _i_ o 1 

i^P2C°p + 27Tixpf + -e-''°u;'Cp^ 



I _|_ 2 1 

—^P2AqP + 2-K'jmxp + -e~"°uj'^Apn 



W„"(0) 
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A'(0) 



\XP 4 J 



1 



1/0/2 



= 'P2 - + 2)P4 - ^P2;^2 - ^AoP2 ) C^P - ^P2CoP ~ ^^XPXP+ 

XP 2 3' 



g (4^Xp)'Ao-poP2 (Co°)2 



+ e 



(47r) 



=''0/2 



+ c^2e--o/2 
:P2i^2 - 



_ XP 2 3 / 

—lp2 - {1 + '2)pa - T,P2V2 - ^AoP2 Aln - lp2Aln - ATrxplm+ 
XP 2 3/ 

1 47r 



Wp{0) 



-finxpAo - nop2 (Aq) 



+ C^2g-.o/2 



= + 7:'^2 + ^Ao ) Apn - Apn 
XP 2 3 ' 



(99) 



IX. APPENDIX II: THE LEAVER SERIES 

Leaver's series is an analytic solution of the outgoing wave of the Regge- Wheeler equation 

f^^[Vi,^ir)~u;']^HW , (100) 
dr* 



where 

'1(1 + 1) ^.2M 



2M 

VBwir) = 1 



r 



(101) 



We will consider the case s = 2, which is for gravitational waves; s = and s = 1 correspond to the massless scalar 
field and electromagnetic waves, respectively. 

Letting g represent the outgoing wave solution, then its analytic expression is 



g[Co, f) = fi+" (f - l)Pe'P'' J2 + ^)nU{s + l + 2p + n,2s + 1, 2pf) , (102) 

11=0 



where r = r/2M and p = —icu = —i2MLU, U is the irregular confluent hypergeometric function, and 

, , T(z + n) , ^ 

is the Pochhammer's symbol. The coefficients a„ are determined by the recurrence relation 

a„a„+i + /3„a„ + 7„a„_i = n=l,2... 

a„ = n < , (104) 

where ao is arbitrary and 

an = (n+ l)(n + 2p + 1) , 

f3n--- [2n^ + {8p + 2)n + 8p^ +4p + l{l + 1) + (1 - s^)] , 

7„ = + 4pn + Ap^ - . (105) 

It must be remarked that, in the derivation of g{uj, f), the time dependence is assumed to be e~"^*, which has a minus 
sign different from our convention e"^*. 

The even-parity perturbation outside the star is described by the Zerilli function, while the Leaver's series gives 
the analytic expression for the outgoing wave solution of Regge- Wheeler equation. Thus, we have to employ transfor- 
mations between the Regge- Wheeler and Zerilli functions in order to apply Leaver's series. The transformations are 
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[/x^ (m^ + 2) + 12iLoM] Z= fj!^ [ii^ + 2) + 



72M2(r-2M)" 
r2 (/iV + 6M) _ 



fliw + 12M 



dr* 



(106) 



[/i^ (^2 ^ 2) - 12iwM] = IJ? {l^^ + 2) + 



72M2(r-2M)" 
r2 (/i2r. + 6M) _ 



Z - 12M— , 
dr* 



(107) 



where /j!^ = {I — + 2) and M is the mass of the star. 
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Modeli^ 


(Tn /ni„ 

^ IL 1 lb 


cr„/m„ 

pi 1 L 




f p 


nnffni 




M/Mp, 


i?(km) 


one 


0.2 


0.5 


2.0 


2.0 


1.3 


0.52 


1.180 


8.657 


two 


0.2 


0.5 


2.5 


2.0 


1.3 


0.741 


1.355 


7.923 


three 


0.2 


0.5 


2.0 


2.0 


2.0 


0.8 


1.178 


7.651 


four 


0.2 


0.5 


2.5 


2.0 


2.0 


1.414 


1.30 


6.841 



Table I. Listed are four specific choices for the parameters of the reiativistic two-polytrope master function given in 
Eq. ([79I). Also listed are the number densities, masses, and radii for the explicit solutions used in the quasi- normal 
mode analysis of the main text. 



Mode 


Re{Auj)/Re{uj) (lO"'^) 


Im{Auj)/Im{uj) (IQ-^) 


1 


3.8 


2.0 


2 


2.1 


2.4 


3 


0.6 


5.4 


4 


1.6 


4.9 


5 


1.8 


9.8 


6 


3.1 


12 



Table II. Listed are the relative changes in the real and imaginary parts of u for the first six w-modes when the step 
sizes are changed from Ar — 10^^ to Ar = 5 x 10^^. 





Model one 


Model two 


Mode 


Re(cjM) 


Im(a;M) 


Rc{ujM) 


Im(cjAf) 


Win 


2.05 X 10-^ 


0.692 


0.151 


0.775 


W112 


0.327 


0.380 


0.428 


0.397 


Wl 


0.497 


0.265 


0.510 


0.194 


W2 


0.848 


0.367 


0.853 


0.316 


W3 


1.183 


0.421 


1.189 


0.369 


W4 


1.514 


0.460 


1.523 


0.406 



Table III. Listed are the frequencies, for I = 2, of the w-modes for models one and two of Table L wm and W112 are 
the strongly damped modes of Fig. 3 with Im(wM) > Re(wM); wi, W2, etc., are likewise the first four modes of Fig. 
3 (reading left to right) whose frequencies satisfy the opposite inequality Im(wM) < Re{ujM). 



Model# 


Mode 


I ^ 2 


I = 3 


I = 4 


one 


f 


0.108 


0.142 


0.166 




Pi 


0.234 


0.281 


0.319 


two 


fo 


0.136 


0.180 


0.212 




fs 


0.157 


0.192 


0.221 




Plo 


0.306 


0.365 


0.414 




Pis 


0.354 


0.417 


0.472 



Table IV. The real parts of the frequencies {ojM) of the f- and pi-modes for / = 2, 3,4 for models one and two of 
Table I. 
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12 3 4 

Re(coM) 

FIG. 3. lm{ijjM) vs Re(ajM) for the w-modes for models one and two described in Table I, setting I = 2. 
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FIG. 4. Ke{mnWn) or Re(mnWp) (since they are indistinguishable) vs r/R for the model one modes wi through W4 listed 
in Table III. 



24 




11 

ll 
1 

1 

1 


1 










1 

1 


1 \\ 
\\ 
\\ 
\\ 
\\ 1 




/ 
// 

r 


\\ 

\ 

N 


/' 


1 

v 










model 1 












one-fluid system 

1,1 



0.2 0.4 0.6 0.8 1 

Re(coM) 

FIG. 6. logj^Q|j4in| vs Re(ajA/), with I = 2, for model one of Table I, and for a one-fluid system composed of neutrons behaving 
as a relativistic polytrope (i.e. /3„ = Pp, cjn = Op, and with the same values for no, /3n, and (t„ as model one). 
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model 2 

one-fluid system 

0.2 0.4 0.6 0.8 1 

Re(coM) 

FIG. 7. logiolAinl vs Re((j;A/), with / — 2, for model two of Table I, and for a one-fluid system composed of neutrons 
behaving as a relativistic polytrope (i.e. /3„ = (3p, an ~ <Jp, and with the same values for no, /3„, and (t„ as model two). 
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FIG. 8. An/n or Ap/p (since they are indistinguishable) vs r/R, with 1 — 2, for the model one f- and pi- modes listed in 
Table IV. 
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FIG. 9. An/n and Ap/p vs r/R, with / — 2, for the model two fo- and fs-modes hsted in Table IV. 
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FIG. 10. An/n and Ap/p vs r/R, with I = 2, for the model two pio- and pis-modes listed in Table IV. 
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FIG. 11. Re(m„Wn) or Re(m„Wp) (since they are indistinguishable) vs r/R, with 1 = 2, for the model one f- and pi-modes 
listed in Table IV. 
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FIG. 13. Re(m„W„) and Re(m„VKp) vs r/R, with / = 2, for the model two pio- and pis-modes listed in Table IV. 
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FIG. 14. Re(//i) vs r/R, with I = 2, for the model one f- and pi-modes listed in Table IV. 
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FIG. 15. Re(ffi) vs r/R, with I = 2, for the model two fs-, pio-, and pis-modes listed in Table IV. 
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FIG. 16. Re(Ar) vs r/R, with / — 2, for the model one f- and pi-modes listed in Table IV. 
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FIG. 17. Re(_R') vs r/R, with Z = 2, for the model two fo-, fs-, Pio-, and pis-modes hsted in Table IV. 
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